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Quantum state reconstruction based on weak continuous measurement has the advantage of being 
fast, accurate, and almost non-perturbative. In this work we present a pedagogical review of the 
protocol proposed by Silberfarb et al, PRL 95 030402 (2005), whereby an ensemble of identically 
prepared systems is collectively probed and controlled in a time-dependent manner so as to create 
an informationally complete continuous measurement record. The measurement history is then 
inverted to determine the state at the initial time through a maximum-likelihood estimate. The 
general formalism is applied to the case of reconstruction of the quantum state encoded in the 
magnetic sub levels of a large-spin alkali atom, ^^^Cs. We detail two different protocols for control. 
Using magnetic interactions and a quadratic ac- Stark shift, we can reconstruct a chosen hyperfine 
manifold F, e.g., the 7-dimensional F = 3 manifold in the electronic-ground state of Cs. We review 
the procedure as implemented in experiments (Smith et a/., PRL 97 180403 (2006)). We extend the 
protocol to the more ambitious case of reconstruction of states in the full 16-dimensional electronic- 
ground subspace {F = 30F = 4), controlled by microwaves and radio-frequency magnetic fields. We 
give detailed derivations of all physical interactions, approximations, numerical methods, and fitting 
procedures, tailored to the realistic experimental setting. For the case of light-shift and magnetic 
control, reconstruction fidelities of ^ 0.95 have been achieved, limited primarily by inhomogeneities 
in the light shift. For the case of microwave/RF-control we simulate fidelity > 0.97, limited primarily 
by signal-to- noise. 



PACS numbers: 32.80. Qk,03.67.-a,03.65.Wj 

I. INTRODUCTION 

An essential tool in quantum information science is 
quantum tomography (QT) [ ]. The ability to estimate 
a quantum state is required to diagnose quantum infor- 
mation processors and evaluate the fidelity of a given pro- 
tocol. The fundamental information-gain/disturbance 
tradeoff in a measurement of a quantum system implies 
that any QT protocol requires multiple, nearly identical 
copies of the state. Typically, this procedure is carried 
out through a series of strong destructive measurements 
of an informationally complete set observables acting on 
repeatedly prepared copies of the system. For this reason, 
QT is generally a time consuming and tedious procedure 
when applied to large dimensional systems [ ] and even 
more so when extended to quantum-process tomography 
in which a whole collection of quantum states must be 
analyzed [3]. 

In some platforms, one has the ability to probe a large 
ensemble of identical systems simultaneously. In this 
case, one can enhance the collection of statistics, e.g., for 
estimation of the probability of occupation in an eigen- 
state under projective measurements. While such projec- 
tive measurements on ensembles can be used for QT [ ], 
in principle one can dramatically improve the speed, ro- 
bustness, and experimental complexity of QT by instead 
employing weak-continuous measurement. In a protocol 
originally developed by Silberfarb et al. [ ], an atomic 



ensemble undergoes a chosen dynamical evolution to gen- 
erate an informationally complete measurement record. 
An algorithm is then used to invert the measurement his- 
tory to determine the maximum likelihood of the initial 
state. For moderately large ensembles, one can attain a 
sufficient signal-to-noise ratio to enable extraction of the 
required information, while simultaneously maintaining 
the quantum projection noise below the intrinsic noise of 
the quantum probe. In this case, quantum backaction is 
negligible and in principle, one can extract the necessary 
data for QT in a single run of the experiment on a single 
ensemble. 

We have previously employed a continuous measure- 
ment protocol to perform QT on the 7-dimensional, 
F = 3 atomic hyperfine spin manifold, in an ensemble 
of cesium atoms [5]. With this tool in hand, we diag- 
nosed the performance of state-to-state quantum maps, 
designed and implemented by optimal control techniques 
[ ]. More recently, this QT protocol was a central com- 
ponent that enabled us to measure the time evolution of 
the quantum state of the spin undergoing the quantum 
chaotic dynamic of a nonlinear kicked top [ ] . Observing 
dynamics of a density matrix for any reasonable duration 
would have been formidable without an efficient method 
for QT at each time step. 

Our goal here is to give a detailed and pedagogical de- 
scription of the protocol, focusing on the ingredients that 
are necessary to enable QT of complex systems, and its 
application to atomic spin systems and other platforms. 
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The remainder of this paper is organized as follows. In 
Sec. II, we give a detailed review of our weak-continuous- 
measurement QT protocol for the general reconstruction 
of a density matrix. In Sec. Ill, we specialize to the case 
of atomic spin systems, and in particular, the control 
of the ground-electronic manifold of cesium atoms. Two 
cases are studied for different control methods: (1) QT 
on the 7-dimensional F = 3 manifold of Cs controlled by 
the quadratic AC-Stark shift and quasi-static magnetic 
fields; (2) QT on the fuh 16-dimensional F = 3 © F = 4 
electronic ground state subspace, controlled by time de- 
pendent radio frequency and microwave magnetic fields. 
In the former case, experiments have been carried out 
and we review those results. The latter case is much 
more ambitious and we outline the challenges as shown 
from our simulations. Finally, Sec. IV includes the con- 
clusions and outlook of this work. 



II. QUANTUM TOMOGRAPHY VIA 
CONTINUOUS MEASUREMENT PROTOCOL 

A. The basic protocol 

The general setting for our protocol is as follows. One 
is given an ensemble of noninteracting, simultane- 
ously prepared systems in an identical state po that can 
be controlled and probed collectively. We seek to find an 
estimate of the state of the system by continuously mea- 
suring some observable Oq. We restrict our attention to 
states in Hilbert spaces of finite dimension d and mea- 
sure traceless-Hermitian observables in the algebra 5u{d) . 
In an idealized form, the probe performs a QND mea- 
surement that couples uniformly to the collective vari- 
able across the ensemble and measures Oq\ For 
a sufficiently strong QND measurement, quantum back- 
action will result in substantial entanglement between 
the particles. For example, such a phenomenon has been 
employed to create spin squeezed states of an ensemble 
when the fluctuations in project ion- valued measurements 
("projection noise") can be seen resolved within funda- 
mental quantum fluctuations in the probe ("shot noise") 
[8-11]. We consider the opposite case of a very weak mea- 
surement such that back-action noise is negligible com- 
pared with the detector noise. In this case, the procedure 
can be analyzed as a single atom control problem in which 
each member of the ensemble evolves under the same dy- 
namics with state p{t). In this case, the measurement 
record is proportional to 

M{t) = TT{Oop{t))^aW{t), (1) 

amplified by the total number of atoms. Here, W{t) is 
Gaussian-random variable with zero mean and variance 
cr^, that is introduced to account for the noise on the 
detector. 

In order to generate a measurement record that can be 
inverted to determine the initial state, one must control 



the dynamics so as to continuously write new informa- 
tion onto the measured observable. To do so, the system 
is manipulated by external fields. The Hamiltonian of 
the system, H{t) = H[(pi{t)]^ is a functional of a set of 
time-dependent control functions, ^i(t), which are chosen 
so that the dynamics produces an informationally com- 
plete measurement record M{t). Since our objective is to 
estimate the initial state of the system from the measure- 
ment record and our knowledge of the system dynamics, 
it is more convenient to carry out the procedure in the 
Heisenberg picture. Expressed this way, the state is fixed 
and control is used to generate new observables that we 
measure. Note, this is generally different from the stan- 
dard Heisenberg picture in that we allow for decoherence 
during the dynamical evolution. We will return to this 
issue below. The measurement record, Eq. (1), is then 
written M{t) = TT{0{t)po)^aW{t). For implementation 
of the algorithm, a time discretization of the problem is 
necessary. We sample the measurement record at discrete 
times so that 

M,=Tt{0^po)^ctW^. (2) 

We have thus reduced the problem of QT to a linear 
stochastic estimation problem. The goal is to determine 
Po given {Mi} for a well chosen {Oi} in the presence of 
noise {Wi}. 

A number of transformations are necessary to increase 
the numerical stability and reliability of the algorithm. 
Let {E^T I /Vd}^ = — 1, be an orthonomal 

Hermit ian basis of matrices, where / is the identity ma- 
trix and Tr(£^Q,) = 0. The unknown initial state, po, can 
thus be decomposed as 

where = Tr(po^Q!) are real numbers. We can then 
write Eq. (2) as 

Mi=Y,^c. Tr((9,^«) + (jWi, or, 

Ol = l ^ ^ 

M = (5r + crW, 

which in general is an overdetermined set of linear equa- 
tions with d^ — 1 unknowns r = (ri, . . . , r^2_i) and where 
= Tr(0,^,). 

Eq. (4) explicitly states that the conditional proba- 
bility of the random variable M given the state r is the 
Gaussian distribution 

V{M.\y) (xexp (^-^(M-(5r)^(M-(5r)^. (5) 

We can use the fact that the argument of the exponent in 
Eq. (5) is a quadratic function of r to write the likelihood 
function 

V{M.\v) cx exp (-^{r - rMLfC-'{r - tml)) , (6) 
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describing a Gaussian function over possible states r cen- 
tered around the most likely state, tml, with covariance 
matrix, C = a'^{0^0)~^ . The unconstrained maximum 
likelihood solution is given by 

TML = ((5^(5)-i(5^M. (7) 

Since we treat the noise as Gaussian, this solution corre- 
sponds to the least squares solution of the linear system 
inEq. (4), [14]. 

The eigenvectors of C ^ specify the directions in the 
operator space $u{d) that have been measured and its 
eigenvalues are the square of the signal-to-noise ratio 
of those measurements. The covariance matrix thus al- 
lows us to quantify the information extracted from the 
measurement record. A sufficient condition for an in- 
formationally complete measurement record (though not 
necessary unit fidelity due to noise) is one for which 
is full rank. This will be true when {Oi} spans 
5u{d). To achieve an informationally complete measure- 
ment record, the quantum system must be controllable in 
the sense that we can map any Oq to any Oi over the Lie 
algebra. At early times during the reconstruction pro- 
cedure, will not be full rank and thus we use the 
Moore- Penrose pseudo inverse in Eq. (7) [ ]. Also, it is 
essential that the dynamics be sufficiently coherent such 
that an informationally complete set of observables can 
be generated before decoherence erases the state. 

B. The positivity constraint 

For an informationally incomplete measurement record 
and/or for finite noise, the unconstrained maximum like- 
lihood solution, Eq. (7), generally produces estimates 
of the density matrix with negative eigenvalues. To ob- 
tain a physical estimate of the state, we must impose the 
constraint that the estimated density matrix be positive 
semidefinite. Such a constraint can be enforced through 
an appropriate parametrization of the unknown initial 
state, e.g., po = T^T where T is a lower diagonal ma- 
trix [1, 4]. Although this parametrization has the advan- 
tage that the estimated state is Hermitian and positive- 
semidefinite by definition, it is not compatible with our 
continuous measurement protocol. A least squares solu- 
tion to Eq. (2) would involve a nonlinear unconstrained 
optimization for which there is no known efficient solu- 
tion. 

We thus turn to constrained numerical optimization to 
find the "closest" positive matrix to the unconstrained- 
maximum-likelihood estimate, i.e., the constrained- 
maximum-likelihood estimate, p. The covariance matrix 
determines a natural cost function metric with which to 
measure the distance these states according to 

llrML - r||' = (rML - rfC-^ruL - r). (8) 

Technically speaking, this quantity is not a norm but 
rather a seminorm only when informationally incomplete 



measurements are considered (C is not full-rank), mean- 
ing that there exist some vectors v such that ||v|| = 
but V 7^ in those cases. The use of this metric can 
be justified as follows. The inverse of the covariance 
matrix, C~^, encodes all of the information about the 
independent directions in operator space that are being 
measured by our procedure. A small eigenvalue of 
means a low signal-to-noise ratio associated with mea- 
surements of the corresponding eigen-operator, and thus 
that little is known about the trace-projection of the ini- 
tial state onto that operator direction. The cost function, 
Eq. (8), takes into account that different directions in 
the space 5u{d) are not measured in the same way and 
weights this in the distance between the initial estimate 
and the positive state. In this way, during the numerical 
optimization, the more uncertain components of p can be 
adjusted more freely than the more certain ones, thereby 
maintaining faithfulness with the measurement record, 
but ensuring positivity. To find the physical estimate we 
thus solve the following optimization problem: 

minimize ||rML - (9) 

subject to 

-/+ ^ raEa>0. (10) 

While there is generally no analytic solution to this prob- 
lem, it takes the form of a standard convex program since 
the matrix is positive semidefinite and both the ob- 
jective and the constraint are convex functions [ ]. The 
optimization is a semidefinite program which is efficiently 
solvable numerically. We implement this in MATLAB 
using freely available convex optimization package, CVX 

[ ]• 

C. Control and dynamics 

An important question that has not yet been addressed 
is the way in which we drive the dynamics to generate 
the measurement record. As discussed above, a sufficient 
condition is that the dynamics generate an information- 
ally complete set of observables {C^i}, meaning that they 
span the Lie algebra 5u{d). The quantum dynamics must 
thus be "controllable" in the Lie algebraic sense [18]. 
The Hamiltonian that governs the dynamics is a func- 
tional of a set of control waveforms such as externally 
applied fields parametrized by frequencies, amplitudes, 
and phases. Our task is then to chose these waveforms 
to generate an informationally complete set of observ- 
ables {Oi} in the desired time. In practice, we fix the 
duration of the measurement record as determined by 
the characteristic time scales for evolution, dictated by 
both the Hamiltonian evolution for the given power in 
the controls and decoherence. We choose the total time 
T to be such that we can attain a good approximation to 
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any unitary evolution matrix in SU{d). The total time 
is then coarse-grained into slices of duration 5t^ consis- 
tent with the slew rates and bandwidth constraints of 
the waveform drivers in the laboratory. We thus reduce 
the problem to specification of a discrete set of waveform 
values compatible with experimental constraints. The 
translation of the discretely sampled parameters to the 
continuous-time waveform depends on the characteristics 
of the physical drivers and the challenges of numerical in- 
tegration, to be discussed below. 

With the Hamiltonian in hand, we must choose the 
control parameters. There is no unique solution; any 
choice that yields an informationally complete set {Oi} in 
the given time series will suffice. In principle, one would 
like to optimize the information gain over time T, and 
yield an unbiased state reconstruction. This amounts to 
optimizing the entropy associated with the eigenvalues of 
the covariance matrix. We have found empirically that 
the landscape for performing such an optimization is not 
favorable, and this approach becomes intractable, even 
for moderately sized Hilbert spaces {d> 9). Instead, our 
numerical studies show that one can achieve the required 
high fidelity and unbiased measurement record by choos- 
ing the control parameters randomly and unbiased over 
a designated interval. We will demonstrate this below 
for the specific example of control and measurement of 
atomic hyperfine spins. A more rigorous justification of 
this approach is still under consideration. We have pre- 
viously seen a connection between evolution via random 
unitary dynamics and the generation of an information- 
ally complete measurement record [ ] . This may give us 
clues to optimally designing the control waveforms. 

An essential component of this protocol is accurate 
modeling of the dynamical evolution of the observables 
measured in the continuous signal. Fundamental to this 
is decoherence induced by the environoment. Under typ- 
ical conditions of Markovian evolution, these dynamics 
are generated by a Lindblad master equation, 

dp{t) 



dt 



= -z[H{t), p{t)] - i ^ {LiL,p{t) + p{t)LlL 



(11) 



The formal solution to this equation is a completely pos- 
itive map on the initial density operator, p{t) = Vt[p(0)] 
with Vt being the solution to 



dV^ 
dt 



(12) 



where T is the time-ordering operator. 

We seek, however, the solution to the Heisenberg evolu- 
tion, given formally by the adjoint map 0{t) = V/[O(0)], 
satisfying 



dt 



(13) 



Naively, one might assume that the generalization of the 
master equation for Eq. (11), to the Heisenberg 

evolution for 0{t) is 



dOjt) _ .t 



dt 



+i[H{t),0{t)] - i ^ {OmlL, + LlL,0{t)) 



Y.LlO{t)L^. 



(14) 



However, this is not generally true since V^oCt ^ CtoVj . 
Note that the correct Heisenberg evolution is 



dO{t) _ , ,t 
dt 



4 [0(0)] 



(15) 



and dO(t)/dt ^ C\[0(t)] unless C is time independent. 
The lack of commutativity between the adjoint map and 
its generator will be the case for a generic time-dependent 
control Hamiltonians under consideration here. Because 
of this, the decohering Heisenberg operators do not sat- 
isfy a time-local differential equation [ ]. This severely 
complicates the efficiency with which we can integrate 
the dynamics to determine the measurement set {Oi}. 

To deal with this problem in moderately large Hilbert 
spaces, as we will discuss in Sec. IIIB, we restrict our 
waveform so that the control parameters are piecewise 
constant over a reasonable duration. Then, over each in- 
terval we can simply exponentiate the Lindbald generator 
of the superoperator map. The traceless operators, O^, 
are "vectorized" to a large column of dimension d^ — 1 
and the superoperator, V^., is a large (d^ ~ 1) x {d'^ ~ 1) 
matrix expanded in the basis E^. Using curved bra- (row) 
ket- (column) notation for the supervectors and superop- 
erators, our integration then takes the form 



{0,\ = {Oo\Vu, 



where 



(16) 



(17) 



For non-piece wise constant controls, this corresponds to 
an Euler integration of the completely positive map. 
Such an approximation will be very inefficient and nu- 
merically unstable for large dimensional systems. For 
this reason a piecewise constant control is best suited to 
our protocol. 



D. Further technical considerations 

Beyond decoherence, an essential ingredient for accu- 
rate modeling of the dynamics is parameter estimation. 
The ability to reach high fidelities for the estimated states 
relies on the assumption that we know exactly how the 
system is evolving at the time that the data is taken so 
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that we know exactly which operators {Oi} are being 
measured. This means that all of the parameters in the 
Hamiltonian must be precisely known before the quan- 
tum state reconstruction is even possible. In practice, 
many such parameters can be precalibrated. However, 
other parameters, such as a background fields may be 
unknown, and other parameters may be inhomogeneous 
across the ensemble. Our protocol is robust because, un- 
like other optimal control tasks, such as state-to-state 
mapping, the exact parameters of the experiment need 
not be fixed. As long as we can determine, post-priori, 
the operating conditions via parameter estimation, and 
the signal is informationally complete, we can extract the 
quantum state with high fidelity. 

An additional parameter we must fix is the initial ob- 
servable being measured. Though abstractly we have 
called this Oq, in practice the true observable may deli- 
cately depend on special alignment of the apparatus. We 
will see how we can use parameter estimation as well to 
fix this observable and the overall calibration of the signal 
in physical units when compared with the dimensionless 
units treated here. 

Finally, one technical detail that we have not discussed 
so far is the signal-to- noise ratio (SNR). Even under ideal 
conditions, the fidelity of the QT is fundamentally lim- 
ited by noise. For Gaussian white noise, it is essential to 
limit the bandwidth in which we analyze the measure- 
ment record and the dynamics must be chosen so that 
the relevant information about the state is contained in 
a narrow frequency band. In addition, 1// noise in the 
detector dictates that the relevant signal be sufficiently 
far from DC. To maximize our SNR for a given ensemble, 
we pass the measurement signal through a narrow band- 
pass filter. In any numerical simulation of the expected 
fidelity, the simulated signal is passed through the equiv- 
alent digital filter as used in the laboratory. This ensures 
equivalent measurement records for equivalent operating 
conditions. 

In the next section, we will apply our QT protocol to 
the reconstruction of states encoded in atomic hyperfine 
spins. After defining the system, we can simulate a mea- 
surement record in the presence of noise, decoherence 
and errors for a given initial state po and use it to run 
it through the algorithm to find the estimate p. In order 
to quantify the performance of our method, we calculate 
the fidelity 



Tr 



'PPo^ 



(18) 



Good performance is judged by high fidelity averaged 
across a collection of randomly sampled states. 



III. APPLICATION TO HYPERFINE SPINS 

In this section, we apply the QT protocol to the recon- 
struction of states encoded in spin systems. Our plat- 




FIG. 1: (Color online) Schematic of our system geome- 
try. A cold gas of atoms is collected from a magneto-optic- 
trap/optical molasses, and optically pumped to form a nearly 
pure ensemble of identical spins. The spins are controlled 
through a combination of light-shift interaction, magnetic 
fields produced by pairs of Helmholtz coils, and microwave 
fields. A measurement of the spins is performed by polariza- 
tion analysis of the transmitted probe. A sketch of the atomic 
level structure for ^^^Cs is shown inset (not to scale). 



form is the hyperfine manifold of magnetic sublevels as- 
sociated with the ground-electronic state of laser-cooled 
alkali- metal atoms, providing a Hilbert space of dimen- 
sion d = {2S ^ 1)(2/ + 1) where 5' = 1/2 is the sin- 
gle valence electron spin and / is the nuclear spin. In 
particular, we work with ^^^Cs, whose nuclear spin is 
/ = 7/2, yielding hyperfine coupled spins of magnitude 
F = 3, 4 and a total Hilbert space of dimension d = 16. 
A schematic of the system, including atomic level struc- 
ture, control, and measurement components, is shown in 
Fig. (1). 

The spins can be controlled through a combination of 
magnetic interactions and off-resonant optical coupling 
from a laser field. The fundamental Hamiltonian is 



H{t) = Al-S- fi-B{t) 



(19) 



where A is the hyperfine coupling constant, fj, is the 
atomic magnetic moment operator, and aij is the atomic 
dynamic polarizability operator, both depending on the 
atomic spin degrees of freedom. Here and throughout, 
we take the laser field complex amplitude, E, to be fixed, 
and control is accomplished through time- variation of the 
magnetic field, B(t). Under typical operating conditions 
where the hyperfine coupling dominates over all other 
forces, the total spin angular momentum, F = /±l/2 = 
3, 4 and its projection along a quantization axis, m, are 
approximate good quantum numbers and define the basis 
of states in the Hilbert space we seek to control. 

A central component of our QT protocol is quan- 
tum controllability. A finite dimensional system with 
a generic Hamiltonian of the form H{t) = Yl^ ^j{t)Hj, 
with external fields determined by Aj(t), is said to be 
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controllable if {Hj} generate the Lie algebra of the rele- 
vant group of unitary matrices on the space [ ]. As we 
generally do not measure the trace of the density oper- 
ator, we restrict our attention to the Lie algebra 5u{d). 
For situations in which we seek control in a single irre- 
ducible manifold F, the relevant algebra is $u{2F + 1), 
F = 3, 4; on the the entire hyperfine manifold the algebra 
is 5u(16). For F > 1/2, this requires Hamiltonians that 
are not linear in all of the components of F. A challenge 
for any implementation is to access nonlinear interactions 
that render the system controllable and coherent. 

Continuous measurement of the system is carried out 
through polarization spectroscopy of a probe laser beam 
that passes through the ensemble while it is being con- 
trolled. The atoms induce a polarization dependent in- 
dex of refraction in a manner depending on their spins' 
state according to the light-shift interaction [ ]. In 
the limit of negligible backaction, the effect of the in- 
teraction is a rotation of the probe's Stokes vector S on 
the Poincare sphere according to the rotation operator 
Ur = exp (-ixo {O) ' S), where xo = ODo{r/2Ac) is the 
characteristic rotation angle depending on the resonant 
optical density, ODq^ and a characteristic detuning from 
resonance, Ac. Taking the z-axis along the direction of 
propagation of the probe, the components of the vector 
of atomic observables that generate the rotations about 
the three axes of the Poincare sphere are. 




0.e3 = ^4V 



F,F' 



A / f f -\- r f 



F'F 



(20a) 
(20b) 
(20c) 



where C^^p are coupling constants that depend on the ir- 
reducible rank-i^ tensor polarizability for the given probe 
detuning /S^^'p from the ground {nSi/2)F to the excited 
{nP'j)F' manifold [ ]. For weak interactions under con- 
sideration here, xo 1, this rotation corresponds to a 
small local displacement. Measurement of the Stokes vec- 
tor component along the direction n then correlates with 
a measurement of the atomic operator n • O. Prepar- 
ing the probe initially linearly polarized along the ei of 
the Poincare sphere, and analyzing along the direction, 
n = cos + singes, the general measurement record 
will be of the form 



M{t) = a {F^Fy 



-h{F,),^cjW(t), (21) 



where a and h are constants that depend on the vec- 
tor and tensor contributions to the polarizability for the 
given detuning, as well as the polarization the analysis 
direction, 0. 

The light-shift interaction also induces dynamics on 
the atomic spin, depending on the probe's polarization. 
The combination of coherent evolution and decoherence 



due to photon scattering can be modeled by a master 
equation of the form [zi] 



dpit) 
dt 



(H,^{t)p{t)-p{t)Hl^{t) 



(22) 



In this equation, projections of operators onto subspaces 
with a given F are denoted, A^^^^ = Pp^APp^ = 
Em,mJ^i''^i)(^i''^il^l^^2,m2)(F2,m2|. The total 
effective Hamiltonian is given by HQf^(t) = Hhf + 
Hsit) + where Hhf is the hyperfine interaction, 

Hb is the magnetic interaction, and the effective (non- 
Hermitian) Hamiltonian accounting for light-shift and 
optical pumping is 



H. 



LS 



eff 



_ 



(€*-DffO(D^'f-^) 



FF' 



^F'F 



-ir/2 



(23) 



Here, Q is the laser Rabi frequency for a unit oscillator 
strength. The strength of the transitions for e-polarized 
light are accounted for by the dimensionless dipole raising 
operator, 

• dJ,,^ = JCj'/ {F'm + q\ Fm; Iq) \F'm + q){Fm\ 

m 

(24) 

where the coefficient JCjp is given in terms of Wigner 
6j symbol 

)Cjp^' = (_i)^W+V(2J' + l)(2F+l)| ""j {-^p] 

(25) 

The Lindblad jump operators are given by 

pi/ a 

describing absorption of a photon with polarization 
€, emission of a photon with polarization and op- 
tical pumping between hyperfine manifolds Fa and 
F5. Transfer of population between sublevels by 
optical pumping occurs at a rate jFama^Fbrnb = 
I {FbTribl Wq^^"" iFaTTia) p. The final term in the mas- 
ter equation, Eq. (22), proportional to p^^^S represents 
transfer of coherences that may exist between hyper- 
fine manifolds, but are preserved in spontaneous emission 
when the detuning of the light is sufficiently large. 

With this general framework in hand, we have the tools 
necessary for our QT protocol: control and continuous 
measurement. We apply this in two examples, consider- 
ing different control scenarios, to achieve quantum state 
reconstruction. 
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A. Reconstructing a single F-manifold via 
light-shift control 

Previous work by [5] and [6] demonstrated the QT pro- 
tocol on the F = 3 ground-state manifold of Cs. In that 
case, the light acted simultaneously as the probe on the 
system and the driver of nonlinear spin dynamics through 
the light shift, and thus made the system controllable on 
SU{2F + 1). For linear polarization of the laser probe 
along X, the effective light-shift Hamitonian Eq. (23) in 
a given manifold F can be expressed in irreducible tensor 
components as [21] 



^(0)_ ^(2)^X^1+1)), ^^^(2)^2 



(27) 

where Ip is the identity operator on the hyperfine mani- 
fold F and 



r/2 



^AF'F + ir/2 



(28) 




FIG. 2: (Color online) Control waveform as a function of time, 
(/)(t), that determines the Zeeman Hamiltonian, Eq. (29). 50 
discrete points of the function were optimized to maximize 
the amount of information acquaried while a cubic spline al- 
gorithm was used to interpolate them. 



are complex coupling coefficients depending on the rank- 
K atomic polarizability. The real part leads to the 
light shift and the imaginary part causes decoherence 
via photon scattering. For emphasis, we have explic- 
itly factored out the characteristic photon scattering rate 
7sc = (l^^r)/(4A^), which sets the time scale for dynam- 
ics on the atom-photon interaction. The real part of p^^^ 
is a shift independent of magnetic sublevel and thus can 
be neglected here. The real part of however, is 

essential for controllability. The best figure of merit for 
nonlinearity vs. decoherence is obtained if the laser is de- 
tuned approximately in between of the excited manifolds 
of the Dl line. Precisely, we choose the detuning from 
(6S'i/2)F = 3 to (6Pi/2)F' = 3 to be Ac/27r = 642.78 
MHz. For this choice of detuning, the nonlinear light 
shift scales relative to ^sc as Re(/3'^^^) = 6.53, while deco- 
herence scales as Im(/3(^)) = -0.23 and Im(/3(2)) = 0.005, 
which implies sufficient coherent control can occur for QT 
before decoherence erases the information. 

We can achieve full control of a hyperfine manifold 
F through a combination of the nonlinear light shift and 
magnetic interaction, as given in the fundamental Hamil- 
tonian, Eq. (19). In the linear Zeeman regime, the 
magnetic moment restricted to this subspace is iip = 
—gpl^B^^ where is the Lande g- factor . Fixing the 
magnitude of B and allowing its direction to vary in the 
x-y plane, the control Hamiltonian is 



HB{t) = QLicos (^{t)F^ ^ sin ^{t)Fy), 



(29) 



where = QfIJ^bB \s the Larmor frequency. The opera- 
tors {Fx^Fy} along with F^, from the light-shift interac- 
tion, form a minimal set of generators of the Lie algebra 
5u(2F + 1) for any F. 

Generation of an informationally complete set of ob- 
servables is achieved through the choice of the waveform 
0(t), the angle of orientation. The Heisenberg dynamics 



evolve the observables according to Eqs. (16) and (17), 
with the Lindbladian given in Eq. (22). As we restrict 
our attention to control of a state in the subspace F = 3 
and the probe is detuned close to the excited state hy- 
perfine splitting in order to generate the nonlinear light 
shift, the population in F = 4 is essentially invisible to 
the probe. Any optical pumping is thus treated as a loss 
and the master equation, Eq. (22), restricted to a single 
manifold, is not trace-preserving. Our reconstruction al- 
gorithm is insensitive to this loss in that we fix the state 
to be normalized through the parametrization, Eq. (3). 
In [ ], the control waveform was designed through a lo- 
cal optimization of the covariance matrix entropy. This 
was possible given the moderate size of the Hilbert space 
{d = 7). For large dimensional system, this approach is 
not tractable, and instead we empirically choose random 
waveforms, as we will discuss in Sec. IIIB. 

As an example of the performance of our quantum 
state reconstruction protocol, we show the reconstruc- 
tion procedure for experimental data in [6]. The fields 
are chosen with nominal parameters Q.l/2it — 17.5 kHz, 
^sc/'^^ = 81-4 Hz. For these time scales, full control- 
lability is achieved in ~ 4 ms. Over this duration, the 
waveform (j)(t) is sampled by 50 control parameters at 
80 /is intervals, consistent with the driver slew rates. The 
optimized angle as a function of time, shown in Fig. 2, is 
made into a continuous waveform via cubic splines. With 
the time-dependent magnetic field interaction, and non- 
Hermitian light-shift Hamiltonian Eqs. (29) and (27), 
respectively, the observables are evolved according to the 
adjoint of the map, Eq. (17). In this case, the integra- 
tion is an Euler approximation to the continuous-time 
differentiable Hamiltonian. Again, because of the mod- 
erate size of the Hilbert space, this is possible. For large 
dimensional systems, such integration is unstable and it 
is essential that the waveforms are piecewise constant. 
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Practical application of the protocol depends on ac- 
curate estimation of all of the parameters that define 
the experiment. Before driving the atoms with the time- 
dependent control waveform that generates the informa- 
tionally complete measurement record, we drive simple 
dynamics so that the Hamiltonian parameters can be 
calibrated. We prepare a spin coherent state, polariz- 
ing the atom along y via optical pumping as our known 
initial state. Fixing the control magnetic field along x, 
continuous measurement of atomic Larmor precession via 
Faraday rotation according the measurement record, Eq. 
(21), with a = 0, is shown Fig. 3. The signal exhibits the 
collapse and revival of Larmor oscillations characteristic 
of the nonlinear spin dynamics first seen in [ 2]. Deco- 
herence via photon scattering is exhibited by the decay 
of the overall signal. In addition, the signal shows in- 
homogeneous broadening due to variations in the probe 
intensity across the ensemble. Finally, technical issues 
such as finite response time affect the time origin of the 
measurement record. All of these features must be accu- 
rately included in our model of the Heisenberg evolution 
in order to obtain high-fidelity QT. 

To perform parameter estimation, we employ a 
least-squares fit between the measured and simulated 
measurement record with cost function C = ||M — 
M(r^L7 7sc7 ^o) IP- Here M is the vector of time-sampled 
data from the calibration run and 

M,(l^L, to) = J {Oi{nL,^7sc,to))f{Od^ (30) 

is the simulated measurement time-series. The unknown 
parameters are the Larmor frequency, I^l, photon scat- 
tering rate, 75c, and origin of time, to. In addition, we ac- 
count for inhomogeneity in the laser intensity through the 
distribution function /(O? where ^ is the ratio between 
the nominal scattering rate and the local scattering rate 
at the position of the atom. The parameter then rep- 
resents the scattering rate at ^ = 1 where /(^) is peaked. 
Note, we take /(^) to be unnormalized. The overall scale 
determines the conversion between the simulated dimen- 
sionless signal and the laboratory measurement record. 

In order to make this inversion tractable, we em- 
ploy a linearly interpolated intensity distribution. We 
parametrize /(^) as a piecewise linear function evaluated 
only at discrete points /(O — + ^n, for ^n-i ^ 
^ < where = [f{^n) - /(^n-i)] / (^n - ^n-i) and 

K = Kn/(^n-l) - ^n-l/(^n)] / (^n - ^n-l)- Making USC 

of this parameterization, the simulated Larmor preces- 
sion signal can be written 



N 



(31) 



where Qn,i = /|;_^ ((9,(^0) and T,,, = 

f^'^ ^ {Oi{^')) d^' . The measurement record can be ex- 
pressed as a linear equation M = Af + aW where 




t[ms] 



FIG. 3: (Color online) Comparison between the experimen- 
tal Faraday rotation signal (blue) and the signal fitted by 
our model (green) for the case of Larmor precession in the 
presence of the nonlinear light-shift. This signal is used to 
calibrate the intensity distribution seen by the atoms. 




FIG. 4: (Color online) Estimated distribution of intensity over 
the atomic ensemble, f{^), where ^ is the ratio of the intensity 
seen by the atoms to the nominal (peak) intensity. 



f = (/07 /17 • • • 7 /at) is the function /(^) evaluated in 
N -\- 1 = 17 discrete points. The least squares solution 
is f = (A-^A)"-*^ A-^M and we use this result as our es- 
timate of the intensity inhomogeneity. This procedure 
is iteratively repeated for different values of nominal in- 
tensity, magnetic field and origin of time until the cost 
function is minimized. 

Figure 3 shows a comparison between the experimen- 
tal signal and the fitted one by the procedure described 
above. The estimated Larmor frequency and peak scat- 
tering rate for this plot are r^i,/27r = 17.469 kHz and 
7sc/27r = 84.7 Hz. Figure 4 shows the piecewise linear 
estimation of the distribution of the intensity of the laser 
probe over the atomic ensemble. We use this distribu- 
tion in the QT run of the experiment to average over the 
intensity inhomogeneity. 
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A final calibration must be performed to determine 
the exact observable we measure in the QT run. For the 
protocol considered here, because we are necessarily de- 
tuned close enough to resonance in order to generate a 
sufficient nonlinear light shift, it is more advantageous 
to measure the birefringence of the atomic sample corre- 
sponding to the observable FxFy-\-FyFx. This observable 
evolves more rapidly into the higher order polynomials in 
F necessary to reconstruct the density operator as com- 
pared to F^, measured in Faraday rotation. Imperfec- 
tions in waveplates, however, imply some uncertainty in 
the exact direction along the Poincare sphere in which 
we perform a polarization analysis. We perform a second 
Larmor precession calibration run, but with a nominal 
quarter waveplate so that the polarization spectroscopy 
corresponds to a birefringence measurement. Accord- 
ing to Eq. (21), the actual measurement operator has 
a contamination of Faraday measurement. Using the 
the estimated distribution of intensity obtained previ- 
ously, /(O? the simulated measured record in this case is 
Mi=aJ {F^Fy + FyF,). f{Od^ + b J (F,), f{Od^. Note 
that in this expression, the unknown parameters are a 
and b. As before, these two parameters are obtained by 
a least squares fit in which we also fit the magnitude 
of the magnetic field and the origin of time. For the 
data shown here a = 0.1613 and b = 0.1598. Though 
the measurement corresponds almost completely to the 
S2 direction on the Poincare sphere, even a small Fara- 
day contamination will produce a noticeable effect in the 
measurement record due to the fact that its measurement 
strength, being a rank-1 tensor effect, is higher than the 
birefringence's, which is a rank-2 tensor effect. Account- 
ing for the small misalignment is essential to guarantee 
high fidelity reconstruction. The result of the simulated 
and measured signals for this calibration step are shown 
in Fig. 5. The excellent agreement between our model 
and the measured signal allows us to achieve high fideli- 
ties in our QT procedure. 

With calibrated parameters in hand, we carry out QT 
of a prepared quantum state. The system evolves under 
the Hamiltonian given in Eq. (29) for an appropriately 
chosen control wave form (/)(t), as discussed above. The 
effect of atomic birefringence on the transmitted probe is 
measured in our polarimeter and the resulting measure- 
ment record is modeled by Eq. (21) with the previously 
estimated parameters. This record is then used in Eq. 
(7) to determine the maximum-likelihood estimate. In 
addition, we carry out a final round of parameter esti- 
mation and fit for the magnetic field amplitudes in the 
two coils, Bx and By^ which produces slightly different 
Larmor frequencies in the x and y directions, and the ori- 
gin of time, to, for the particular run. While the B-field 
is unlikely to change between the initial calibration and 
the QT run, we find that refitting is necessary to ensure 
high fidelity. The origin of this discrepancy is unknown, 
but might be explained as a compensation for incorrect 
fitting of the intensity distribution and its complicated 
effect on the overall dynamics. Removing the need to fit 
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FIG. 5: (Color online) Comparison between the experimental 
birefringence signal (blue) and the signal fitted by our model 
(green) for the same Larmor precession dynamics given in Fig. 
3. This signal is used to calibrate the measurement basis. 



for intensity inhomogeneity would greatly simplify the 
QT protocol, as we consider in future generations of our 
protocol, discussed in the next section. 

Figure 6 shows the fitted and the experimental recon- 
struction signal of a simple example: a spin coherent 
state along y. Once this fit is done, Eqs. (9) and (10) are 
used to find the closest positive-semidefinite estimate of 
the initial state. Figure 7 shows a bar plot of the density 
matrix elements of the initial state and those of the recon- 
structed one. Qualitatively, we see how similar these two 
plots are and quantitatively we calculate the fidelity, Eq. 
(18), to be about 0.95. To compare this with the average 
expected behavior, including much more complex states 
produced in optimal control, we simulate measurement 
records based on our master equation, added noises, and 
signal processing. For 1000 initial states randomly sam- 
pled using the Hilbert- Schmidt measure, we estimate our 
reconstruction fidelity to be about 0.998, limited only by 
the noise in the probe. We attribute the lower fidelity 
of the experimental example above as arising from limi- 
tations in our ability to determine (estimate) all the pa- 
rameters present in the actual system dynamics. Relax- 
ing some of these limitations will increase the fidelities we 
can reach in experimental situations. In the next section 
we describe our next-generation protocol that addresses 
some of these limitations while simultaneously allowing 
for reconstruction of the full 16-dimensional Hilbert space 
associated with the electronic ground-states. Realization 
of this protocol would represent a substantial advance in 
complex quantum state reconstruction. 



B. Reconstructing the entire hyperfine- manifold 
via rf and microwave control 

In this section, we develop the tools necessary for re- 
constructing the entire 16-dimensional hyperfine ground 
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FIG. 6: (Color online) Comparison between the actual recon- 
struction signal (blue) and fitted by our procedure (green) for 
an initially prepared spin coherent state along y. 



(a) 



(b) 




FIG. 7: (Color online) Comparison between the actual initial 
state (a) and the state reconstructed by our procedure (b). 
The absolute values of each element of the density matrix are 
plotted as a function of its index. The fidelity achieved by 
our procedure in this case was ^ 0.95. 



state manifold (F = 3 F = 4) of cesium. To achieve 
this we make use of the type of control developed in [ ] , 
which employs microwave (/iw) and radio-frequency (RF) 
modulated external magnetic fields to drive the atoms. 
The RF-fields drive rotations on the F = 3 and F = 4 
manifolds, whereas the microwave fields drive a reso- 
nant transition between two Zeeman levels in F = 3 
and F = 4, thus making the system fully controllable. 
Since the light-shift interaction is not necessary for con- 
trollability, in principle, the control can be achieved in 
a decoherence free way [ ]. However, in practice, since 
a laser is used to measure the system, some decoherence 
will be present that must include in our model. 

The fundamental Hamiltonian, Eq. (19), can be writ- 
ten in this case H{t) = Hq -\- Hrf + H^u) + Hls- The 
free Hamiltonian, Hq, includes the hyperfine interaction 
and the Zeeman shift produced by the bias magnetic field 
interaction, which is necessary to define the quantization 



axis and to resolve microwave- induced transitions of a 
pair of magnetic sublevels. Keeping terms to quadratic 
order in the field strength. 



X 

2 V* ' ^ 



(P4-P3) + l^o(Fi4^-5.Fi3)) 

(32) 



where ujhf/'^^ ^ 9.19 GHz is the the hyperfine split- 
ting for cesium, Pp, is the projector operator onto the 
F = 3,4 manifolds, VLq = gAjUBBo is the Larmor fre- 
quency produced by the bias field Bq in the z direc- 
tion, X = {ge — gi)/^BBo/ujHF with gj ^ —0.0004 and 
ge ^ 2.0023 being the nuclear and electronic g- factors, 
respectively [zej]. The Lande g- factors for each manifold 
F have opposite signs and are given by gs ^ 0.2499 and 
g4 ~ —0.2507 and its ratio gr = Igs/g^l ^ 1.0032 is the 
relative g-factor between the F = 3 and F = 4 manifolds. 
The last term in Eq. (32) is the quadratic Zeeman shift, 
where a = x'^ujhf/{'^I + 1)^- 

The RF-control Hamiltonian, Hrf^ is produced by a 
RF-magnetic field that oscillates in the x and y directions 



Hrf = n,{t) cos {ujRFt - Mt)){F^^^ 
^ny{t) COS {ujRFt-Mt)){F^''^ - 



(33) 



where the Larmor frequency is defined as ^i{t) = 
g4/J^BBi{t), for i = X,?/, (I)x{t), and (l)y{t) are the con- 
trol phases of the magnetic field in the x and y directions 
respectively and ujrf is the frequency at which the RF 
fields are modulated. This Hamiltonian allows for inde- 
pendent SU(2) rotations within the F = 3 and F = 4 
manifolds. 

For control via application of microwave radiation, we 
consider a purely cr+-polarized field. While in prac- 
tice the polarization of microwaves is not well controlled 
at the position of the atoms, this is not critical, since 
ultimately we will drive a selected two-level transition 
through its unique resonance frequency. Nonetheless, we 
allow for the effect of off-resonant ac- Stark shifts caused 
by microwaves and choose one polarization to analyze, for 
simplicity. Under this assumption, the microwave control 
Hamiltonian is 



Hfj^w =^fiw{t) cos {ujf^wt - (l)fiw{t))x 

3 



^ (4,m + l|3,m;l,l)4-), 



(34) 



m= — 3 



where 

^flW (^) 



the bare microwave Rabi frequency is 
= IJ^sB^y^it)^ (l)fiw{t) is its control phase. 



(4, m + 1| 3, m; 1, 1) is the Clebsch-Gordan coefficient 
associated to the transition |3,m) |4, m + 1), 

and cri"^^ = |4,m + 1) (3,m| + |3,m) (4,m + 1|. This 
Hamiltonian couples Zeeman levels in the two different 
manifolds, taking into account the resonant as well as 
the off-resonant transitions. 
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With the control Hamiltonian in hand, we perform em- 
ploy the rotating wave approximation (RWA) to elimi- 
nate the explicit time dependence. Following Appendix 
A, the free Hamiltonian, Eq (32), written in the rotating 
frame and after the RWA is given by Eq. (A2) 



where cr^^^ = |3, m) (3, m| — |4, m ■ 
and (jy are given in Appendix A. 



1) (4, m + 1|, and (Jx 



ffo = (^(l 



25 

?r) + ya ) {Pa 



Ps) 



+ no{l-gr)Fi^^-a{{Fi^yf-{F('^f), 



(35) 



where we consider the resonant case Arf = A^^ = 0. 
The RWA for the microwave Hamiltonian is straightfor- 
ward since in general Q/j^^j <C ojhf- However, we must 
pay special attention to keep the correct off-resonant 
terms that lead to microwave-induced AC Stark shifts. 
The resulting microwave Hamiltonian in the RWA is 



^flW (^) 

2 

"/iW ( 

SUJRF 



[cos (j)^^ {t)cFx + sin (j)^^ ifycfy] 
^Iwit) ^ |(3,m;l,l|4,m + l)p 



-crl 



m 



The case of the RF-Hamiltonian is much more com- 
plex. In order to maintain rapid control, the RF-Larmor 
rotation frequencies must be sufficiently large. However, 
if the bias field is not sufficiently large, the condition 0.^^ 
<^ OJRF will not be fulfilled. In this case, we must 
consider higher order corrections to the RF Hamiltonian 
in the RWA. To do this, we follow [ ], and employ the 
method of averages. The details of this procedure are 
discussed in Appendix A. The resulting RF-Hamiltonian, 
Eq. (A23), is 



HRF{t) = 



2 
2 

E 



cos 



cos0y(t) (fW 



(1 + 9r)9r 



i?(3) 



■sin (fW 



(3 - gr)9r 



(37) 



(l-2cos2</.i(t))(FW 



— sm 

OUJRF 



frFi'^: 



where, as before, we consider the resonant case and set 
ujRF = ^0- This differs from the Hamiltonian given in 
[23] in two ways. We account for the relative magnitudes 
of the ^-factors in the upper and lower manifolds due to 
the small nuclear magneton, ^ 1. Additionally we 
maintain terms of order Q'^/ujrf^ ^ = x^y^ which lead 
to Bloch-Seigert-like shifts and extra corrections due to 
counter-rotating terms. 

We must also transform the light-shift Hamiltonian, 
Eq. (27), to the rotating frame. The RWA is an excel- 
lent approximation in this case since generally we will 
have parameters such that the Zeeman splitting is much 
larger than the rate of coherent coupling induced by the 
light shift, fto ^ Re(/3^^^)75c. Making the appropriate 
unitary transformation and averaging over the rapidly 
varying terms, as described in Appendix A, we obtain 
the effective light-shift Hamiltonian in each manifold F, 



lie 



fLS 



eff,F 



I] 



6 



(38) 

Note, under the RWA, the light-shift does not drive co- 
herences between magnetic sublevels defined by the quan- 



tization axis. Such coherent couplings are no longer reso- 
nant in the presence of a strong bias field. Also note that 
in considering dynamics over the full hyperfine manifold, 
we must retain the real part of the scalar light shift, since 
generally ^ /3^^ . The scalar contribution to the light 



shift thus drives coherences between the F 
manifolds. 



3 and F = 4 



We now have all the ingredients to proceed with our 
simulations. The system evolves according to the mas- 
ter equation, Eq. (22), with the effective Hamiltonian 
expressed in the RWA and the control Hamiltonian de- 
scribed above. In addition, special care must be taken 
when considering the full master equation. All operators, 
including the Lindblad jump operators, must be written 
in the rotating frame and the RWA should be applied ac- 
cordingly. This is essential in order to account for spon- 
taneous emission processes that become distinguishable 
once the energy degeneracy is broken by the shift pro- 
duced by bias field. The RWA in the master equation 
is achieved by explicitly calculating the transformation 
W {t)Wq^^'^U{t), and averaging the superoperator map 
over the rapid oscillations. 

Following [ ], full controllability of the system can 
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FIG. 8: (Color online) Randomly sampled control waveforms 
that determine the phases of the applied RF and microwave 
magnetic fields as given in Eqs. (36,37). These waveforms 
produce high fidelity estimates with reasonable stability. 



be achieved by keeping the RF-Larmor and /iw-Rabi fre- 
quencies, and constant in time, while vary- 
ing the control phases (l)x{t), ^yif)i ^^^(t). Due to 
the size of the Hilbert space, finding a set of control wave- 
forms is a very challenging task. Optimizing the entropy 
of the Gaussian probability distribution, as mentioned in 
Sec. Ill A, is generally an intractable problem. Instead, 
we choose the control waveforms as piecewise random 
functions. We have found in our numerics that this is 
sufficient to generate an informationally complete mea- 
surement record. For the RF (/iw) waveforms, the phase 
is chosen uniformly between — tt and tt and kept constant 
over intervals of 30/is (20/is). After a total time of 2 ms 
we ensure that we have sufficient information in the mea- 
surement for QT. Fig. 8 shows an example of the control 
phases used in this work. In general, almost all the con- 
trol waveforms designed in this random way will produce 
informationally complete measurement records for most 
initial states. However, numerical stability during the 
use of Eq. (7) becomes an issue for certain waveforms. 
We thus choose a set of waveforms that produce the most 
stable results by repeating the design procedure several 
times. 

To complete our protocol, we must choose the detuning 
Ac of the probe. In the previous section, this detuning 
was chosen to maximize the nonlinear light shift relative 
to photon scattering. In the current context, we have 
much more flexibility, since full control of the Hilbert 
space can be achieved without the light shift. There are, 
however, many technical considerations that inform the 
choice of detuning. Firstly, the measurement strength is 
proportional to 75c, so we can never make the measure- 
ment record free of decoherence by detuning further off 



resonance [^ 3]. Moreover, at very large detunings, in or- 
der to maintain a reasonably large 7^0, we would need a 
large probe intensity for which shot-noise-limited detec- 
tion is difficult. For these reasons, the light-shift-driven 
dynamics must be included in the analysis. An additional 
technical issue is the effect of inhomogeneity in the light 
intensity across the ensemble. Indeed, the difficulty in es- 
timating the distribution of intensities caused substantial 
complexity in the reconstruction algorithm, as discussed 
in Sec. HI A, and ultimately limited the fidelity of the 
protocol. Mitigating this effect would greatly improve 
the performance. 

From Eq. (27), we see that the scalar, Re(/3'^^^) and 
tensor, Re(/3^^^), light shift components are responsible 
for the introduction of inhomogeneity in the problem. 
While there is no choice of detuning that makes both 
terms exactly zero, the scalar light shift will cause the 
largest problems, and our goal is to cancel it. Putting all 
of these considerations together, we choose a relatively 
small detuning where the measurement strength can be 
large at low intensity. For such a detuning, only one 
F manifold is effectively coupled to the light, and the 
other is so far from resonance that its coupling is neg- 
ligible. We chose this to be F = 3, and detune of the 
Dl line, (651/2)^ = 3 to (6^1/2)^' = 3. Detuned within 
the hyperfine splitting of the excited state, we can find 
a "magic wavelength" at which Yle{l3^^) = 0. Using the 
tensor coefficients given in [ U], we find the magic detun- 
ing on the Dl line is Ac/27r = 291.89 MHz. The resid- 
ual light shift is due to the tensor term, proportional to 
Re(/33 ) = 1.35, which together with microwave and RF 
fields, drives the spin dynamics during the course of the 
measurement. 

To study the performance of our protocol, we per- 
formed numerical simulations of the expected measure- 
ment signal for a random pure state sampled from the 
Haar measure. We choose the following parameters: 
nx/27T = ny/27r = 15 kHz, n^^/27T = 33 kHz, Qo/^tt = 
1.0 MHz, and a characteristic photon scattering rate 
75c/27r = 410.7 Hz. Furthermore, we add Gaussian white 
noise to the signal so that the signal-to-noise ratio is 
100. As a first check, confirm that our choice of detun- 
ing makes our system more robust to light-shift inhomo- 
geneity, by adding Gaussian fluctuations in the intensity 
across the ensemble, and then averaging the result. For 
example, if we choose the same detuning as used in Sec. 
Ill A, Ac/27r = 642.78 MHz, f3^^^ 7^ 0, and we see that the 
average signal leads to fidelities < 0.80, whereas a similar 
simulation with the optimized detuning, Ac/27r = 291.89 
MHz produces fidelities > 0.90. Fig. 9 shows qualita- 
tively a comparison between the simulated measurement 
records, averaged over a Gaussian distribution of inten- 
sities for these two detunings and the simulated signals 
with a fixed, nominal value of intensity. It is clear that 
the optimized detuning produces much better results, 
making both averaged and nominal signals look very sim- 
ilar. This will simplify our procedure for estimating the 
intensity distribution seen by the atomic ensemble. A fit 
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FIG. 9: (Color online) Simulated measurement signal for a 
random pure state chosen from the Haar measure for differ- 
ent choices of detuning of the laser probe. (Top) Ac/27r = 
291.89MHz, the "magic wavelength" at which the scalar light 
shift is set to zero for F = 3. The red line is an averaged sig- 
nal over a Gaussian distribution of intensity. The green line is 
the signal that a system would produce if it evolves under the 
nominal value of intensity. (Bottom) Ac/27r = 642.78MHz, 
the detuning used in previous experiments. The red line is 
an averaged signal over a Gaussian distribution of intensity. 
The green line is the signal that a system would produce if 
it evolves under the nominal value of intensity. The signal is 
more robust at the "magic wavelength". 



to a Gaussian distribution will be sufficient to capture 
the effects of the inhomogeneous light shift. 

A number of calibration errors are possible in this sys- 
tem. We anticipate the need of fitting similar parameters 
to the ones in the last section, i.e., the RF and /iw mag- 
netic field amplitudes as well as the origin of time. Least 
squares techniques, similar as the ones discussed in de- 
tail in Sec. Ill A, should suffice to achieve good accuracy 
parameter estimation and high-fidelity reconstructions. 
Finally, we choose to measure Faraday rotation instead 
of birefringence of the probe, which further simplifies the 
protocol since we no longer must fit for the initial mea- 
surement operator. Our simulations below show that this 
choice leads to high fidelity QT. Thus, for the simulations 
shown in this section, we have chosen our initial observ- 
able to be the Faraday operator (9o = • 63 as given in 
Eq. (20c). 

We have run several simulations to test the perfor- 
mance and efficiency of this protocol. We numeri- 
cally generated a measurement record for different initial 
states and different noise realizations and amplitudes ac- 
cording to Eq. (2). Then a Bessel bandpass ffiter from 
6 to 80 kHz was applied to the simulated measurement 
record in order to limit the noise in frequency compo- 
nents that are not present in the measurement. The 





FIG. 10: (Color online) Real and imaginary parts of the 
elements of po, the initial state used in the simulations. 
This is a non-trivial state used for illustration, \ip) = 

^ V^ig'*^ + V^cat^^ /a/2 5 consisting of an equal superposition 

= 



of a spin squeezed state in the F = 4 manifold 

exp {— iO.SF^^} \F = 4:, nix — 4) and a "cat state' 

(|F = 3,m;, = 3) + |F = 3,m^ = -3))/V2, in the F = 3 
manifold. 



same filter is applied the Heisenberg picture observable 
0{t) to account for all dynamical effects the signal un- 
dergoes. Once this is done, Eq. (7) is used to find the 
unconstrained maximum likelihood estimate of the initial 
state and the convex program Eq. (9) is solved to find 
the physical density matrix that best represent the mea- 
sured data. In order to quantify the performance of our 
method, we calculate the the fidelity between the initial 
and estimated state, Eq. (18). 

As an example, we simulate the reconstruction 



position of 
manifold. 



))/V2, 

consisting of an equal super- 



( 



^cat 



the nontrivial state, 

shown in Fig. 10, 

spin squeezed state in the F = 4 

exp {-z0.5F|} |F = 4,ma, = 4) and 
a "cat state" in the F = 3 manifold 



€1) 



(|F = 3,m^ = 3) + |F = 3,m^ = -3))/V2. In Fig. 11, 
we show how our procedure converges as a function of 
time to an estimate of the initial state with high fidelity. 
Within the few first microseconds of the simulation there 
is little information, and the protocol returns the maxi- 
mally mixed state as the estimation of po- However, as 
the time passes, more information about the information- 
ally complete set observable is acquired and the protocol 
makes better guesses of the initial state. For a SNR of 
100 simulated here, within 1 ms, a fidelity of > 0.97 is 
achieved. 

In Fig. 12, we show the fidelity between the estimated 
state and the initial state for a random pure state in the 
16 dimensional Hilbert space as a function of time for 
different signal-to-noise ratios. Although this plot shows 
the performance of the reconstruction protocol for a par- 
ticular state, the same behavior is seen for most random 
states sampled from the Haar measure. In fact, for 200 
of such random pure states, we achieved an average fi- 
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Re{p) 



Im{p) 



(a) 



t = 0.06ms 
J=' = 0.0625 




FIG. 11: (Color online) Real and imaginary parts of the esti- 
mated initial state p. (a) After 60/is not enough information 
has been collected, thus there is no convergence of our algo- 
rithm, and the maximally mixed state is guessed, (b) and 
(c) More information is acquired as time passes and higher fi- 
delities are obtained, (d) A high fidelity estimate is obtained 
after 1.5ms of simulation. Slightly higher fidelities are seen 
for longer simulation times. 



delity ~ 0.977 with standard deviation of 
2ms and an SNR of 100. 



0.006 after 



IV. SUMMARY AND OUTLOOK 

In this work we have presented a comprehensive re- 
view of a protocol to perform fast, robust, high-fidelity 
quantum tomography (QT) based on continuous mea- 
surement of an informationally complete set of observ- 
ables. This procedure is applicable when one has access 
to a large ensemble of identically prepared systems in a 
product state, collectively coupled to a probe field. For 




- SNR=30 

- SNR=60 

- SNR=100 



0.2 0.4 0.6 0.8 1 

t[ms] 



FIG. 12: (Color online) Fidelity, Eq. (18), of the recon- 
structed state respect to the actual, random, pure state as a 
function of time. A jump in the fidelity is seen at t ~ 0.8ms, 
before an informationally complete measurement set is ob- 
tained, due to the positivity constraint. Shown is the per- 
formance of the QT protocol under different signal-to-noise 
ratios (SNR). A fidelity of - 0.977 is seen for a SNR of 100. 



weak measurement back-action, the probability distribu- 
tion of parameters that define the density matrix, con- 
ditioned on the measurement record, is Gaussian, and 
the problem maps onto one of classical stochastic pa- 
rameter estimation. With sufficient signal-to-noise, the 
density matrix can be found from a single measurement 
record using a maximum-likelihood estimate. A physical 
density matrix, with positive eigenvalues, is then found 
using a convex optimization algorithm that searches for 
the closest positive state to the unconstrained maximum- 
likelihood prediction. 

We have applied this protocol to the problem of QT of 
hyperfine spins in cold atomic ensembles. A key com- 
ponent of our procedure is to drive the system with 
well-chosen control fields so as to generate an informa- 
tionally complete set of observables over the course of 
the measurement record. We have presented here two 
methods for achieving this, one based on combinations 
of magnetic-generated Larmor precession and a nonlin- 
ear spin rotation generated by the light-shift, and an- 
other approach based on combinations of microwave and 
radio- frequency-driven spin rotations. The former allows 
us to reconstruct the density matrix associated with one 
hyperfine manifold F, as has been demonstrated in ex- 
periments on the 7-dimensional F = 3 manifold of Cs 
atoms [vj]. The latter is more ambitious, allowing us to 
reconstruct the entire electronic ground state subspace 
(16 Hilbert space dimensions for Cs). 

Our protocols rests on the assumption that once the 
measurement record is obtained, we can invert its his- 
tory to estimate the initial quantum state of the ensem- 
ble. It is thus essential to accurately model the atomic 
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dynamics and measurement of the observables, includ- 
ing known sources of imperfections. We have given a 
detailed discussion of the master equation governing the 
atomic dynamics and a measurement model based on po- 
larization spectroscopy. With these in hand, we showed 
how to simulate the measurement record, with empha- 
sis on limitations, challenges, and the steps needed to 
make our protocol reliable and stable. For the case of 
light-shift control, we tested our state-of-the-art proto- 
col with experimental data and found that we could effi- 
ciently reconstruct the 7-dimensional quantum state with 
a fidelity 0.95, limited primarily by the difficulty in ac- 
counting for the effects of the inhomogeneous intensity. 
Removing this constraint, our simulations show that we 
should obtain a fidelity >0.99 with a SNR of 100. In 
the case of full control on the 16-dim spin system via RF 
and microwave driving, we simulated noisy measurement 
records and used these as inputs to our reconstruction al- 
gorithm. We found that we can rapidly achieve average 
fidelity >0.97 at the same SNR for a measurement time of 
only 2 ms. These initial results bode well for high-fidelity 
reconstruction of a quantum state in such a large Hilbert 
space. Experimental studies are underway. With such a 
tool, we can explore the implementation of qudit unitary 
transformations for quantum information processing [24] , 
and nontrivial dynamics as seen in quantum chaos ['^J. 

Our approach to quantum tomography could be im- 
proved in a number of ways. While we have found that 
random waveforms are sufficient for generating informa- 
tionally complete measurement records, the nature of op- 
timal waveforms (in time and/or average fidelity) remains 
open. Additionally, as the protocol has analogies with 
classical stochastic estimation, we see potential for im- 
proving the reliability and stability of the reconstruction 
procedure by employing the data processing tools such 



as Kalman filters [27] and other methods of estimation 
theory, which we plan to explore in future studies. 
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Appendix A: RWA corrections 

This appendix describes in detail the derivation of the 
Hamiltonian that governs the dynamics hyperfine mani- 
fold driven by RF-and microwave magnetic fields as well 
as light shift interactions. In general, for applicable pa- 
rameters, we must go beyond the usual linear Zeeman ef- 
fect and first order rotating wave approximation, which 
substantially complicates the Hamiltonian beyond that 
presented in [ ] . To begin we transform to a frame that 
is rotating at the frequency of the control fields, accord- 
ing to the unitary transformation U (t) = UrfU^uj^ where 

Urf = exp [-iconFt{F^^^ - ))] (Ala) 
Of 

U^w=exv[-i-{P^-Ps)] (Alb) 

with 6 — uo^w — (^4 + 1^3)^ RF^ where and ms label 
two Zeeman levels corresponding to the F = 4 and F = 3 
manifolds respectively. It then follows from Eq. (32), 



U\t)HoU{t) 



dU 
~dt 



-(1 



, 25 
9r) + ytt 



{IArf - A^^) ) (P4 - P3) - A^^Fi^) + {Arf + l^o(l - 9r)) Fi') 



.a((Fi^))^-(Ff))^), 
(A2) 



where we have chosen 7714 = 4, ms = 3, A^^^? = ujrf — ^q-, 
^fiw (^i^w - ^0, with cjo (^HF + (4 + 3^r)^o + 7a 
being the on-resonant transition for the two-level system 
formed by the stretched states |3, 3) and |4, 4). Although 
our goal is to be as close to resonance as possible {Arf = 
=0)5 in practice we must also account for nonzero 



detunings that might result, e.g., from gradients in VLq 
across the ensemble. 

Going to the rotating frame, the RF Hamiltonian in 
Eq. (33) transforms to H'j^p{t) = U\t)HRF{t)U{t), 
yielding 
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H'nAt) = ^(cos {2ujRFt - Mt)) + cos {Umi^'^ - 9rF^^^) - ^(^^^ ^^^RF^ - M^)) + sin ((/>.(t)))(F^(^) + gr 

+ ^(cos {2uRFt - Ut)) + COS ((/>,(t)))(F^(^) - ^.Ff )) + ^(sin (2c^^^t - ^y(t)) + sin ((/>,(t)))(Fi^) + gr 

(A3) 

I 

In the same manner, we can write the microwave Hamil- 
tonian, Eq. (34), in the rotating frame H'^^(t) 



-(cos (2cJ^^t - (l)i^w{t)) + COS (0^^(t)))cr^ + 



(sin (2cj^^t - (l)^w{t)) + sin (^^w(t)))cr^ 



^/iw (^) 



^ (3, m; 1, 1| 4, m + 1) (cos (2cj^^t + 2(m - 3)cJi?Ft - (|)^,w{t)) + cos (2(m - 3)ujRFt + ^^^(t)))cr^^^ 



7717^^3 



^/iW (^) 



^ (3, m; 1, 1| 4, m + 1) (sin {2uj^^t + 2(m - 3)wfiFt - <P^,w{t)) + sin (2(m - 3)ujRFt + 4>^.w{t))Wy'^ 



(A4) 



r 



where a 



(m) 



-i|3,m) (4,m + l| + z |4, m + 1) (3, m|. 
Note that we have explicitly separated the resonant terms 
from the off-resonant ones. The off-resonant interaction 
produces a AC-Zeeman shift of the magnetic levels that 
must be accounted for in the regime we consider in our 
simulations. We have defined for the resonant transition 
c^x = and cFx = cFy . Finally, the effective light-shift 
Hamiltonian, given by Eq. (23) written in its irreducible 
tensor representation (as in Eq. (27)), can be expressed 
in the rotating frame as 



rrfLS 
^eff 



(Fi^) cos (tORFt) + F^^^ sin {ujRFt)f 



(A5) 



where is given in Eqs. (28). 

Given the Hamiltonian in the rotating frame, we pro- 
ceed to apply the RWA. Typically this is a straightfor- 
ward task, equivalent to dropping the rapidly oscillating 
counter-rotating terms. For the case RF Hamiltonian, 
since the Zeeman splitting induced by the basis mag- 
netic field, 1^0 may not be much larger that the driving 
frequency, urf, a second order correction of the RWA 
is needed. In the remainder of this appendix, we follow 
[26] and use the method of averages for ordinary differ- 
ential equations, which we briefly review, to provide the 
required correction to the RWA. 

Given a set of first order differential equations of the 
form ^ = e f(x, t,e), where x represents the state of 



the system, f(x, t,e) is a periodic function with period 
T, and e is a small parameter, we seek an approximate 
solution of the equivalent averaged system y under the 
transformation x = y + e a;(y,t,e), where a;(y,t,e) is 
also periodic with period T. The averaging theorem says 
that the equations of motion of the equivalent system 
are dy/dt = e f(y) + fi(y,t,e) + 0{e^), where f(y) = 
^/o ^{y^t,e)dt 



Noting that when only the RF part of the Hamiltonian 
is present in the problem, the complete dynamics of the 
system can be described in the SU(2) group, and thus, 
all the dynamics of the system can be described by the 
Heisenberg equations of motion of F^^ Fy and F^. We 
carry out this calculation for the F = A and F = 3 man- 
ifolds separately since there is no coupling between them 
in the absence of the microwaves. Moreover, we assume 
a small enough bias field Bq so that we can neglect the 
quadratic Zeeman shift introduced in Eq. (32); for a very 
large bias the standard RWA is sufficient. For illustra- 
tion, we discuss in detail the second order correction to 
the F = 4 manifold RF Hamiltonian. 



In this case, it is convenient to define the small pa- 



rameter e = cq/cjOrf where eo = y l^^(t) + ^l{t) to allow 

the RF Larmor frequencies to be different. Turning the 
microwave Hamiltonian off {fl/j^w = 0) and neglecting the 
second order Zeeman shift, the Hamiltonian restricted to 
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the F = 4 manifold can be written 



H^'\t) = (^(1 - 9r) + ^(7A^^ - A,^)) 



(A6) 



where 



x{t) 



+ 



^^(COS (2LjRFt - Mt)) + COS iMt))) 

".w... ^^^^ 



-(sin(2cJi?F^ - + sin (^^(t))), 



and 



v{t) 



+ 



^^(sin {2ajRFt - Mt)) + sin (Mt))) 



(cos {2ujRFt - 4>yit)) + COS i4>yit)))- 



(A8) 



The only terms that participate in the averaging process 
are the fast oscihating ones, while the slow varying terms 
are treated as constant. 

We now proceed to calculate the function uj{y,t') = 

Jq (f (y, t") —i{y))dt" by again integrating only over the 
fast varying terms 



(A14) 



where 



X(i) 



(sin {2L0RFt - (j)x{t)) + sin {(j)x{t))) 



+ 



( - COS {2cORFt - 4>y {t) ) + COS {ct>y{t))), 



eo 



(A15) 



The Heisenberg equations of motion for the compo- 
nents of the total angular momentum can then be easily 
written 



(4) 



dF, 



(4) 



dt' 



2 ^ eo 



dF, 



(4) 



dt' 



= e 



2 

x{t') 



( 

eo / 



t;(t') 



_ 

2 ^ 2 



(A9a) 
(A9b) 
(A9c) 



where we have scaled the time so that t' = ojRFt. This 
system of differential equations is in the form needed to 
apply the averaging theorem when we note that 



, f(x,t') 



v{t')F 



(4) 



2AFi") 
-2AFi^) 



-x{t')F^'^ 
x{t')F^^^ - vi^F^'^ 

(AlO) 

where, for convenience, we have defined A = Arf/cq. 
Transforming the original system to the averaged equiv- 
alent one, we have 



and 



+ 



^^(cos {2ujRFt - Mt)) - cos(0^(t))) 



-(sin {2ujRFt - (t)y{t)) + sin {(t)y{t))). 



We can thus write 



A 



fi(y,^0^ 





A 




+ 4 








XF 



(4) 



TF, 



(4) 



(A17) 



where 



A 



1 m) 



cos {24)x{t)) 



cos {2,j>y{t)) 



(A18) 





r F^^) 1 




" uFi") + 2AFi') ■ 




y ^ 


Fi^) 


-xFi^) - 2AFi^) 


(All) 













where 

X ^ 

and 



^- cos((/),(t)) + ^^sin((/)^(t)). 



eo 



^-sin((/),(t)) + ^^cos((/)^(t)). 



eo 



eo 



(A12) 



(A13) 



eo 



^cos(0,(t)), 
eo 



(A19) 



and 



^^^^)-cos(0.(t)) + ^sin(0,(t)). 



eo 



eo 



(A20) 



Putting all together, the Heisenberg equations for the 
components of the total angular momentum, up to second 
order correction of the RWA, are 
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d_ 



X 










e 

~ 2 










XH'^ - vF^'' 





A 



y 

(4) 



A 



(4) 



(4) 
z 



(A21) 



(4) 



Equivalent ly, using Eq. (A21), we can write the Hamil- 
tonian, Eq. (A6), up to second order correction 



+ (^cos(</>.(t)) + ^sin(</>.(t))) +^ (^^^^ 

^ (^sin(</>.(t)) + |^cos(</>.(t))) + ^^ (cosiMt)) ^^'^^(Mt)))) 

+ " ^ (2^- W)) + W (l - 2 cos (2(/)^(t))) + 2Q,{t)ny{t) sin - F^^\ 

(A22) 



Using a similar procedure to the one detailed above, a 
second order correction for the Hamiltonian acting on 
the F = 3 manifold can also be obtained. Putting all the 



second order correction terms together, the RF control 
Hamiltonian for the full ground manifold, in the RWA, 
corrected up to second order can be written 



HRF{t) 



2 

2 

2 

2 
1 

16w_Ri; 

9r 
16u)RF 



cos(</.,(t))(FW 
cos(,^,(t))(FW_5^ 



f^o(l - Qr) 



(sin {Mt))Fi'^ - 9r cos )) 



^ _ l^o(l-gr) ^ ^(3) 

2corf 



sin {(j)x{t))F^ ^ 



— sm 
Arf 



y I I sin(0^(t)) 



2UJRF 

(sin(</.,(t))Ff 



grSm{(j)y{t))F^ ^ 

lit) (l - 2 cos (2(/),(t))) + (^-y _ 2 cog (2(/)^(t))) + 2n,{t)ny{t) sin - ^y{t)) ] F^^^ 



1 - 2 cos ( 
\ - 2 cos ( 



1 — 2 cos 
1 — 2 cos 



lit) (l - 2 cos {2Mt))) + W (l - 2 cos (2(/)^(t))) - 2n,{t)ny{t) sin - j F^ 



(A23) 



r 
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